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Abstract. A method based on the simple shear flow modified by the introduction of a deterministic non-conservative force 
and a stochastic process is proposed to measure the Navier-Stokes shear viscosity in a granular fluid described by the Enskog 
equation. The method is implemented in DSMC simulations for a wide range of values of dissipation and density. It is observed 
that, after a certain transient period, the system reaches a hydrodynamic stage which tends to the Navier-Stokes regime for 
long times. The results are compared with theoretical predictions obtained from the Chapman-Enskog method in the leading 
Sonine approximation, showing quite a good agreement, even for strong dissipation. 



1. INTRODUCTION 

Under the assumption of molecular chaos, the velocity distribution function of a dense fluid of inelastic hard spheres 
obeys the Enskog kinetic equation, suitably modified to account for the collisional inelasticity through a constant 
coefficient of normal restitution a | 1]. By an extension of the Chapman-Enskog method, the Navier-Stokes transport 
coefficients for a single fluid have been derived in the first Sonine approximation as nonlinear functions of the packing 
fraction <j) and the coefficient of restitution 01 . The shear viscosity for a dense binary mixture has also been evaluated 
within the same Sonine approximation 0. 

Given that the first Sonine approximation is known to be rather accurate in the case of elastic hard spheres |4], a 
natural question is whether this approximation is still reliable for finite inelasticity. This point is of interest since in 
the Sonine approximation the distribution function is expanded around the local equilibrium distribution. However, 
a granular system is never in equilibrium |5]. In the undriven case, the uniform reference state corresponds to the 
so-called homogeneous cooling state (HCS), in which the time dependence only occurs through the temperature. The 
velocity distribution function of the HCS deviates from a Gaussian, as measured by a non-zero value of the fourth- 
cumulant and by an overpopulated high-energy tail |6]. This fact could cast some doubts about the accuracy of the 
Sonine approximation to compute the transport coefficients. 

The shear viscosity T] is perhaps the most widely studied transport coefficient in a granular fluid. For a low-density 
gas (0 — > 0), Brey et al. |7] measured rj from DSMC by analyzing the time decay of a weak transverse shear wave, 
and observed a qualitative good agreement with the Sonine prediction. However, while this method is physically quite 
elegant, it is not perhaps sufficiently efficient to measure tj accurately. More recently, the low-density Navier-Stokes 
transport coefficients of the HCS have been measured from DSMC by using the associated Green-Kubo formulas 
|0]. In the case of a system heated by the action of an external driving or thermostat, the associated shear viscosity 
(which is different from that of the HCS) has been computed from the Chapman-Enskog method in the first Sonine 
approximation and from DSMC, for a dilute single gas |9] as well as for a dense binary mixture On the other 

hand, to the best of our knowledge, the shear viscosity of the HCS for a dense gas described by the Enskog equation 
has not been measured from DSMC before. This issue is the main goal of this paper. As a second point, we want to 
address the question of whether a granular fluid reaches a hydrodynamic regime that can be described by the Navier- 
Stokes equations, even far from the quasi-elastic limit. In fact, a granular fluid in an inhomogeneous steady state is 
inherently non-Newtonian la. Hill , so that a Newtonian regime requires unsteady states sufficiently separated from the 
initial preparation. 

In order to investigate the problem, we have carried out numerical simulations of the Enskog equation by means of 
an extension of the DSMC method 1 12] for 0.6 < a < 1 and < <f> < 0.5 and have considered the so-called simple shear 



flow state. To allow the granular fluid to approach a Newtonian regime, an external "friction" force with a negative 
friction coefficient is applied to compensate for the inelastic cooling, so that viscous heating produces a monotonic 
decrease of the Knudsen number (reduced shear rate) with time. In addition, to mimic the conditions appearing in 
the Chapman-Enskog method to Navier-Stokes order, a stochastic process is introduced, according to which every 
particle has a certain probability per unit time to have its velocity replaced by a new velocity sampled from the HCS 
distribution. 

Our simulation results confirm that, regardless of its initial preparation, the fluid reaches after a few collisions per 
particle a hydrodynamic state whose temporal evolution is governed by that of the granular temperature. Moreover, 
when the shear viscosity is measured as the ratio between the shear stress and the shear rate for long times, it agrees 
reasonably well with the theoretical expressions derived in the first Sonine approximation \2], the deviations being 
comparable to or smaller than the ones in the elastic case. 



2. MODIFIED SIMPLE SHEAR FLOW 

The simple shear flow is macroscopically characterized by a constant density «, a uniform granular temperature T, and 
a linear velocity field u(r) = a • r, where the rate of strain tensor a is ay = a5 !V 5, y , a being the constant shear rate 1 13]. 
At the level of the one-body distribution function, the simple shear flow becomes uniform in the Lagrangian frame of 
reference, i.e., /(r, v,f ) = /(V(r),f ), where V(r) = v — u(r). Consequently, the Enskog equation in this problem reads 

d t f-aV yl —f+F[f]=J[f,f]. (1) 
oV x 

Here J[f,f) is the Enskog operator, which in the simple shear flow is given by 

J[fJ] =o 2 X{n) JdVi | fl fa©(&-g)(a-g)[a- 2 /(V')/(V' 1 +a.CT)-/(V)/(V 1 -a-(7)], (2) 

where a is the diameter of a sphere, %{n) is the contact value of the pair distribution function, O is a unit vector directed 
along the centers of the two colliding spheres at contact, g = V — Vi is the relative velocity, a < 1 is the coefficient 
of normal restitution, a = (JG, and the primes denote pre-collisional velocities, namely V' = V-j(l + o; _l )(ff-g)ff 
and Vj = V i + j(l + a _1 )(c • g)G. Finally, F[f] in Eq. is an operator representing some type of external action 
on /, absent in the true simple shear flow problem, that will be chosen later on. This extra term is assumed to verify 
the conditions 

J d\{l,\,V 2 }F[f] ={0,0,-3nyT/m}, (3) 

so that it does not affect the mass and momentum conservation laws, but in general contributes to the energy balance 
equation through the "heating" rate y. The time evolution for the granular temperature is 

d t T = -^ a p xy -;T+yr, (4) 

in 

where P xy — + PjL is the total shear stress and £ is the cooling rate, 

Ifj = J dVmViVjf(V), P}j = l -±^mo 3 % J dadidjj d\ J dVie(ff-g)(d.g) 2 /(V + a ■ ff)/(Vi), (5) 

C = ^(l-« 2 )z| £ /ff|«'v| £ /V 1 ©(a-g)(a-g)V(V + a.C7)/(V 1 ). (6) 

The first term on the right-hand side of Eq. ® represents viscous heating effects, the second term corresponds to the 
cooling due to the inelasticity of collisions, while the third term is the contribution due to the external action F[f]. 

In the true simple shear flow, i.e., with F[f] — 0, the temperature evolution is governed by the competition between 
the viscous heating and the collisional cooling effects until a steady state is reached when both effects cancel each 
other. However, this steady state is inherently non-Newtonian full , so that the Navier-Stokes shear viscosity coefficient 
cannot be measured, even in the quasi-elastic limit. The domain of validity of the Newtonian description is restricted 
to small values of the Knudsen number Kn = X/l, where A = 1 / '\f2Kna 1 '%{n) is the mean free path and I is the 



characteristic hydrodynamic length. In the present problem, the only hydrodynamic length is I = y/2T /m/a. In the 
steady shear flow Kn is a function of a that cannot be controlled independently. In order to reach asymptotically small 
values of Kn for any value of a we need to avoid a steady state and have a monotonically increasing temperature. 
Consequently, we must modify the original shear flow problem by introducing an external driving mechanism, 
represented by F[f], which exactly compensates for the collisional cooling term in Eq. 10}. Specifically, the heating 
rate introduced in Eq. (0 is chosen as 7 = £, so d t T = — (2/3n)aP fy . 

So far, apart from the requirement 7 = £, the explicit form of F[f] remains still open. Here we will determine 
F[f) by requiring that the kinetic equation Q to first order in Kn be the same as the one found by applying the 
Chapman-Enskog method for states close to the (local) HCS |2]. To that end we formally expand / as 



/(V)=/(°)(V)+/«(V)- 



(7) 



where p k > (V) is of order k in Kn. Moreover, the time dependence of / occurs entirely through the temperature. Note 
that, by definition, a ~ Kn. On physical grounds, Pxy is at least of first order in Kn. Therefore, given that 7 = £, 
d t T ~ Kn 2 and so d t f can be neglected to first order. Inserting Eq. |7} into Eq. {0, we get, through first order in Kn, 



(0)1 



: y(OWO)j(0)i 



(8) 



where 



F\fW] ~aV y — /(°) = -L\fW]+aA y df^/dV x , (9) 
' dV x ■ L J 

J^[X,Y] = o 2 %(n) jdVi J da®{a-g){a-g)[a- 2 X{\')Y{\' y )-X{\)Y{\ x )], (10) 

L[X] = -/%(°),X]-/(°)[X,/(°)], (11) 



MX] = a 3 Z («) / rfVj / £ /a©(a-g)(a-g)a i [a- 2 / (0) (V')^(V' 1 )+/ (0) (V)X(V 1 )] 



We require that / (0 ) be the HCS. This implies that | 



2 s «?V 



V/(°)(V) 



(12) 



(13) 



where £(°) is the cooling rate in the HCS, which is obtained from Eq. (|6j by setting a = and / — * f(°\ Next, we 
consider the linear integral equation (|9jl for /W. This equation coincides with the one derived in Ref. from the 
general Chapman-Enskog method specialized to Vn — VT = V • u = 0, provided that 



F[f^] = -C {0) Td T f 



(i) 



I/-(o)A 
2 s dV 



v/ ll, (v) +k (0) / (1) (v) 



(14) 



where in the last step we have taken into account that Kn °c T x l 2 and, by dimensional analysis, /^'(V) = 
«Vq 3 4>(V/vo)Kn, where vq = \J2T /m is the thermal velocity, <I> being a dimensionless function of the reduced ve- 



locity. The simplest choice for F\f] compatible with Eq. (0 (with 7 = Q and with Eqs. Jl 3i and dl4i is 



F[f\ 



(15) 



In summary, our modified simple shear flow problem is described by Eq. Q along with the choice d!5l >. The first 
term on the right-hand side of Eq. (I15> represents the effect of a deterministic non-conservative force of the form 
j£mV, which does work to compensate for the collisional energy loss. The shear viscosity of a granular fluid mixture 
when only this term is present in F[f] has been determined from the Chapman-Enskog method and measured in 
DSMC simulations ylLUJ]. However, this term does not suffice to get the Navier-Stokes shear viscosity of the HCS, 
but it must be supplemented by the second term on the right-hand side of Eq. J15I . The latter term represents the action 
of a stochastic process, whereby each particle has a probability per unit time equal to j £ of changing its velocity by 
a new velocity sampled from the (local) HCS distribution /(°) . When this stochastic term is moved to the right-hand 
side of the Enskog equation Q, it can be interpreted as a BGK-like relaxation term. 
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FIGURE 1. Left panel: Parametric plot of the reduced normal stresses PfjnT versus Kn 2 for a = 0.8 and = 0.5. Right panel: 
Parametric plot of the fourth cumulant c, relative to that of the HCS, versus Kn 2 for a = 0.9 and = 0.2, 0.3, and 0.4. 



3. DSMC RESULTS 

We have solved the Enskog equation Q with the presence of F\f] as given by Eq. (1151 by means of an extension to the 
Enskog equation II 211 of Bird's DSMC method. Since the system is uniform in the Lagrangian frame, there is no need 
to split it into cells. In the implementation of F[f] the cooling rate is self-consistently measured in the simulations. 
The stochastic term is simulated by randomly choosing during a small time step 8t a fraction \ C,dt of particles; the 
original velocity V id of each one of those particles is replaced by a new velocity V new = VTC, where C is the velocity 
of a particle in a reservoir kept at the HCS normalized to unit temperature. 

The shear rate and the initial temperature are such that Kn = 0. 1 at t = for all density and coefficient of restitution. 
Besides, the initial velocity distribution is a Maxwellian. In the course of the simulations, the kinetic and collisional 
transfer contributions l|5} to the pressure tensor are evaluated. From the shear stress P xy (t) the shear viscosity is 
measured as a function of time as 77(f) = —Pxy(t)/a. Since the Knudsen number Kn(f) « 1/ y/T{f) monotonically 
decreases with increasing time, the Navier-Stokes shear viscosity is identified as 77(f) for long times. 

As said in Section^ our main objective is to compare the kinetic (r\ k ) and total (rj) shear viscosity measured in the 
simulations with the expressions derived from the Chapman-Enskog method by using the first Sonine approximation 
/W — > —a{mr\ k /nT 2 )fMV x V y , /m being the Maxwellian distribution. The expressions are [2] 

*, 1 \ l-§0y(0)(l + «)(l-3«) 
77 K (a,0) = tj -| — — — — — — — - - — , (16) 
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(17) 



where 770 = (5 / '16a 2 ) y/ 'mT / '% is the shear viscosity of a dilute gas in the elastic limit, (j) — §«<7 3 is the packing 
fraction, and co is the fourth cumulant of the HCS. A good estimate of cq is |6] 

. , 32(l-a)(l-2a 2 ) 
C0(a) = 81-17a + 30a 2 (l-ar (18) 

In addition, we take the Carnahan-Starling approximation %{$>) = (1 — 0/2)/(l — (j)) 3 . 

Figure[2shows the kinetic part of the normal stresses P^(t), relative to nT(t), and the fourth cumulant c(f ), relative to 
its HCS value cq, as functions of the time-dependent Knudsen number Kn(f ) = A /£(T) = a(T/l2<j>x(<l>)y/T(t)/m, for 
some representative cases. Note that, as time grows, the Knudsen number Kn monotonically decreases from its initial 
value Kn = 0.1, behaving as Kn(f) ~ f _1 for asymptotically long times. We observe that in the long-time limit, i.e., 
for Kn — > 0, the system tends to an isotropic state with a fourth cumulant c equal to that of the HCS. This supports the 
expectation that the asymptotic state of our modified simple shear flow is the HCS, despite the fact that the temperature 
is increasing rather than decreasing, in agreement with Eq. It is worth remarking that, according to Fig.^ after a 
short transient period the fluid reaches a hydrodynamic regime where the normal stresses and the cumulant are linear 
functions of Kn 2 (Burnett-order effects). 
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FIGURE 2. Parametric plot of the reduced kinetic shear viscosity rj k (a,(f>)/ri k (l,(j>) versus Kn 2 for a = 0.8 and <j> = 0.2. The 
dotted line represents the estimated Navier-Stokes value. 



As an illustration of how the Navier-Stokes shear viscosity is evaluated from DSMC, Fig. |2 shows the time- 
dependent kinetic shear viscosity r\ k {t) = —P k (t)/a, relative to its (time-dependent) theoretical Sonine value in 
the elastic limit, as a function of Kn 2 (f) for the case a = 0.8, (j) = 0.2. Figure [2] clearly shows that the ratio 
ri k (a,(j))/r] k (l,<j)) reaches a plateau for long times (small Knudsen numbers) that can be identified as the Navier- 
Stokes value. The same procedure has been followed to measure the kinetic and collisional parts of the Navier-Stokes 
shear viscosity for different values of dissipation and density. 




FIGURE 3. Plots of the kinetic part of the shear viscosity (left panel) and of the total shear viscosity (right panel), relative to their 
respective theoretical values in the elastic limit in the first Sonine approximation, as functions of the packing fraction for a = 0.9, 
0.8, and 0.6. The symbols are simulation results, while the lines are the theoretical predictions in the first Sonine approximation. 




FIGURE 4. Plots of the kinetic part of the shear viscosity (left panel) and of the total shear viscosity (right panel), relative to 
their respective theoretical values in the elastic limit in the first Sonine approximation, as functions of the coefficient of restitution 
for 0=0, 0.2, and 0.4. The symbols are simulation results, while the lines are the theoretical predictions in the first Sonine 
approximation. 



The density dependence of r\ k and 77 for three values of the coefficient of restitution is displayed in Fig. [3] while 
Fig-Sshows the influence of dissipation for three values of the packing fraction. We observe a general good agreement 
between DSMC results and the predictions from the first Sonine approximation, even for strong dissipation and high 
densities. Both theory and simulation show that, at a given value of a, r] k (a,(j)) > r\ k {\,ty) if the packing fraction 
is smaller than a certain value 0q(oc), while r] k (a,<j)) < r] (l,<j>) if <j> > <j>rj(<x). A similar behavior occurs for the 
total shear viscosity with a different value <j>o{<x). The influence of a on both 0q and 0q 1S rather weak; according 
to Eqs. (OB-©, O ) = (0.12,0.09), (0.13,0.09), and (0.15,0.10) for a = 0.9, 0.8, and 0.6, respectively. As 
a consequence, while in a dilute granular gas (0 < 0.1) the shear viscosity increases with inelasticity, the opposite 
happens for sufficiently dense fluids (0 > 0.1). 



4. CONCLUDING REMARKS 

In this paper we have proposed a method to measure the Navier-Stokes shear viscosity of a moderately dense granular 
fluid described by the Enskog equation. The idea is to consider the simple shear flow (which is uniform in the 
Lagrangian frame), modified by the presence of a deterministic non-conservative force (which compensates for the 
collisional cooling) along with a stochastic BGK-like term. Under these conditions the Knudsen number of the problem 
decreases with increasing time, so that the system reaches a hydrodynamic Navier-Stokes regime for long times. This 
procedure allows one to evaluate the Navier-Stokes shear viscosity in an efficient way by means of the DSMC method. 
The simulation results have been compared with predictions from the Chapman-Enskog expansion in the first Sonine 
approximation |2]. The results show that the Sonine predictions compare quite well with the simulation data for the 
wide range of dissipation (a > 0.6) and density (0 < 0.5) explored. This agreement is significant, given that, in contrast 
to the elastic case, the reference state is not a (local) Maxwellian but the (local) HCS,which exhibits non-Maxwellian 
features, such as a non-zero fourth cumulant and an overpopulated high-velocity tail 1 6] . It is interesting to remark that 
the accuracy of the Sonine approximation found here is comparable to the one observed for elastic fluids, so that one 
could expect to improve the agreement by considering more terms in the Sonine polynomial expansion. Finally, we 
want to stress that the system, after a short transient period, achieves a hydrodynamic stage prior to the Navier-Stokes 
regime. This stage could be used to study nonlinear transport properties (e.g., shear thinning and viscometric effects), 
although this issue is beyond the scope of this paper. 
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